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Abstract 

A new approach known as flat histogram method is used to study 
the lb J Ising spin glass in two dimensions. Temperature dependence 
of the energy, the entropy, and other physical quantities can be easily 
calculated and we give the results for the zero-temperature limit. For 
the ground-state energy and entropy of an infinite system size, we 
estimate e° = -1.4007 ± 0.0085 and s° = 0.0709 ± 0.006, respectively. 
Both of them agree well with previous calculations. The time to find 
the ground-states as well as the tunneling times of the algorithm are 
also reported and compared with other methods. 
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1 Introduction 



The equilibrium properties of spin glass have remained a great challenge in 
numerical simulations. Investigating the equilibrium ground-state structure 
of spin glass is also important and interesting. In the last 20 years, there has 
been a great deal of work on spin glass It is generally agreed that the 
simplest spin glass system for most theoretical work is the Edwards- Anderson 
(EA) model, whose Hamiltonian is 

H= - J^3^^^j^ (1) 

<«J> 

where cTj takes on the values ±1 and the sum goes over the nearest neighbors. 
The Jij are dimensionless variables which describe the random interactions 
between the spins and are taken as Jij = ±1. In two dimensions, a phase 
transition occurs only at zero temperature 0, ^, Q for this kind of ± J Ising 
spin glass with nearest neighbor interactions. This model has been stud- 
ied previously by the transfer matrix method 0, ^, replica Monte Carlo 
method 0, |^ , multicanonical ensemble method [0 and many other methods 
(see ref [0 for a review). 

The traditional Monte Carlo methods mostly concentrate on generating 
standard statistical ensembles, e.g., the canonical ensemble or microcanonical 
ensemble. Using the canonical ensemble simulations, we need to simulate at 
different temperatures to get full information about the system. It is tedious 
to calculate certain thermodynamic quantities like the free energy and the 
entropy since the density of states cannot be obtained directly from the simu- 
lation data. The correlation between subsequent configurations generated by 
canonical ensemble simulations also causes the ergodicity problem for some 
systems. In 1991, Berg proposed the multicanonical ensemble method ^ 
to overcome the above shortcomings of simulations on canonical ensemble. 
The multicanonical ensemble is an ensemble where the probability P{E) of 
having energy E at equilibrium is a constant. The multicanonical method 
has been very successful in solving the systems that involve energy barriers. 

Recently, Wang proposed a dynamics which can generate a flat his- 
togram in the energy space as the multicanonical method. This dynamics 



has some connections with the broad histogram method [|1^, which does not 
give the correct microcanonical average 0]. Similar to the broad histogram 
method, the new dynamics is also based on {N{a, AE)) , the (microcanon- 
ical) average number of potential moves which increase the energy by AE 
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in a single spin flip. A cumulative average (over Monte Carlo steps) can be 
used as a first approximation to the exact microcanonical average in the flip 
rate. Thermodynamic quantities can be then calculated from the simulation 
data with ease. In this paper, we use the new method to study the ther- 
modynamics as well as ground-state properties for the two-dimensional Ising 
spin glass system. 

In Section 2, the flat histogram transition matrix Monte Carlo dynamics 
is described. Using the flat histogram sampling, we get the average number 
of potential moves {N{a, AE)) e, which can be used to construct a transi- 



tion matrix Monte Carlo dynamics in the energy space |Tl|. We apply the 



new method to two-dimensional Ising spin glass and present some numerical 
results in Section 3. In the last section, we give a conclusion to the new 
method. 



The transition matrix Monte Carlo dynam- 
ics with the flat histogram sampling 



To connect our dynamics with single-spin-flip Glauber dynamics [|T2|, we 
restrict the protocol of each move to be single-spin flip in the following dis- 
cussion. For a given state a with energy E, consider all possible single-spin 
flips. The single-spin flips change the current state into possible new 
states, with new energy E' = E + AE. For two-dimensional Ising spin glass, 
AE = 0, ±4, and ±8. We classify the new states according to AE and 
count the number of N{a, AE). Since each move from the state a of energy 
E to the state a' of energy E' and the reverse move are both allowed, the 
total number of moves from all the states with energy E to E' is the same 



as from E' to E. Thus, we have ITS 



E N{a,AE)= Y: N{a',-AE). (2) 

E{a)=E E{a')=E+AE 

The microcanonical average of a quantity A{a) is defined as 

(3) 

where the summation is over all the configurations having energy E and 
n{E) is the density of states. In terms of the microcanonical averages, we 
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can rewrite Eq. as 



n{E){N{a,^E))E = n{E + ^E){N{a' ,-^E))e+ae. (4) 

Eq. is the basic result of the broad histogram method |]TU[. While the 
broad histogram random walk algorithm is not correct, Eq. (P is not prob- 
lematic and taken as the starting point of the flat histogram sampling. 

We select a site to flip at random. The flip rate for a single-spin flip from 
state a with energy E to a' with energy E' = E + AE is chosen as 

Then the detailed balance condition for this rate 

r{E'\E)P{a) = r{E\E')P{a') (6) 



is satisfied for P(o") oc l/n(E(a)). Thus the energy histogram is fiat |13 



PiE)= PMcxn(E)^ = const. (7) 

Since {N{a, AE)) e is not known in general, an approximation scheme 
should be used to start the simulation. For those E which we have not visited 
yet, we simply set r{E'\E) = 1. Then a cumulative average (over Monte Carlo 
steps) can be used as an approximation to the exact microcanonical average 
in the flip rate. We have numerical evidence that this procedure converges 
to the exact result. 

We can then construct a transition matrix Monte Carlo dynamics in the 
energy space [Tl| with {N{cr, AE)) e- For a single-spin-fiip Glauber dynamics 
with energy change AE, the flip rate is given as 

1 A F 

u>(AB) = -[l-tanh(— )]. (8) 



Since there are (on average) {N{a, AE))e different ways of going from E to 
E' = E + AE, the total probability for transition from E to E' is 

W{E + AE\E) = w{AE){N{a,AE))E, forAE^O. (9) 
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The diagonal elements can be determined by J2ae W{E + AE\E) = 1, since 
the total probability from E to E' is 1. This new dynamics in the space of 
energy E is related to single-spin-fiip dynamics by |Tl[] 

WiE'\E) = -±- J: J: r(aV). (10) 

""^-^^ E{u)=E E(u')=E' 

where T{a'\a) is the transition matrix of the single-spin-fiip dynamics. The 
equilibrium state of the transition matrix gives the canonical probability 
distribution of energy Pt{E) oc n{E) exp^—E/ksT). 

An important aspect of this dynamics is that we can calculate the thermo- 
dynamic quantities easily by just performing one simulation for each coupling 
state Jij. The density of states n{E) can be obtained through Eq. (H). Once 
we have the density of states n{E), we can obtain Pt{E) and then calcu- 
late any thermodynamic quantities of interest. In actual implementation, we 
usually determine Pt{E) directly from the detailed balance equation 

W{E + AE\E)Pt{E) = W{E\E + AE)Pt{E + AE) (11) 

instead of solving Eq. (^. From Eq. (P), we know, the transition matrix 
W{E'\E) can be formed at any temperature once the quantity {N{a, AE))e 
is computed accurately. In other words, the Monte Carlo computation is 
uncorrelated to thermodynamics. The temperature dependence enters only 
after simulation in the weighting formula. 

Like Berg's multicanonical ensemble simulations, our dynamics also gen- 
erate a multicanonical ensemble in the energy space. From this point, both 
of the two dynamics have the same goal of flattening the space of energy. But 
they are quite different in implementation. In the multicanonical ensemble 
method, the flip rate is chosen as the inverse of the density of states n{E), 
parametrized in some way. To start the simulation, we need give an estimate 
of n{E), since n{E) is not initially known. Thus, the efficiency of this method 
is determined by the goodness of the estimated n{E). If n{E) is not given 
properly, say far off the true density, the simulations may get stuck in some 
region. With our method, we sample the energy space with a flip rate which 
is related to the density of states through Eq. (^). The central quantity is 
{N{a, AE)) E which can be quite accurate in a short simulation time. And 
the accuracy of this quantity is improved by further simulations. We then 
provide an alternate for the estimate of n{E), which leads to a more efficient 
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way for simulating the multicanonical ensemble. 



The flat histogram also generalizes easily to multi-variate models . An 
example is the Ising Spin Glass model with overlap parameter q which has 
the Hamiltonian 

^ = - E JrJ^H - E JrJ^'^' - h E • (12) 

<i,j> <hj> 



E refers to the first two interaction terms involving the coupling constants Jij. 
In this bivariate case, the quantity {N{a, AE)) generalize to {N{a, AE, Aq)). 
It can be easily shown that the detailed balance condition is now 

n{E, q){N{a, AE, Aq))E.g = n{E+AE, q+Aq){N{a', -AE, -Aq)) E+AE,g+Aq 

(13) 

with n{E,q) as the new "density of states". The algorithm gives a flat 
histogram in both E and q. 



3 Numerical results 

We have performed simulations on lattices of size L = 4, 10, 16, 24 and 32. 
Each simulation starts with independent random numbers. To illustrate the 
performance of our algorithm, we define the time as the average number 
(over coupling constant Jij) of Monte Carlo steps needed to reach the ground- 
states. A Monte Carlo step is defined as flipping each spin on the lattice once 
(on the average). Table 1 gives an overview of typical time in Monte Carlo 
steps to reach the ground-states, starting from an arbitrary energy level. The 
time to reach the ground-states depends on the size of the system and also 
the random interactions. We consider a large number of random coupling 
states to make the statistical error small enough in Table 1. The simulations 
are long enough to ensure that the ground states are really reached. In Fig.l 
we plot the time tl versus lattice size L on a double log scale. The data are 
consistent with a straight-line fit, which gives the finite-size behavior 

Tl oc L^-'^\ MC steps. (14) 

The corresponding CPU time for a Digital Alpha 600M workstation is also 
shown in Table 1. For accuracy, 5 independent runs are performed for each 
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lattice size to obtain the average CPU time. Up to L = 32 the CPU time 
can be approximated by a polynomial function of L^ ^*^. 



L TLjMC Steps) CPU time(second) 

4 45.33 ± 0.91 0.003 ± 0.001 

10 2752 ± 127 0.43 ± 0.012 

16 25761 ± 1504 8.55 ± 0.55 

24 216884 ± 18444 189 ± 15 

32 759609 ± 80310 733 ± 86 



Table 1: Average MC steps and CPU time to find the ground-states for 
different lattice size. 



We also consider the tunneling time which is defined as the average Monte 
Carlo steps needed to move from Emax to Emin, or from Emm to Emax- Note 
that Emin is the same as ground state energy and Emax is N — Emm- We 
note that Berg's definition about tunneling time is slightly different from 
ours. During the simulation, Berg imposed a constraint Jij = 0. But for 
our method, both the time and the tunneling time will not be affected 
significantly by the imposition of the constraint. We start the simulations 
from an arbitrary energy level. Table 2 gives an overview of the tunneling 
time obtained using the two methods. The power law fits are 

tm.c. oc and tf.h. oc L^•°^ (15) 

for Berg's method and current fiat histogram method, respectively. It shows 
that they basically give the same tunneling time. 



We also compared with Hatano's result |15| that autocorrelation time 
scales approximately as volume N of the system. We look at the tunneling 
time which is a better measure of the algorithm's efficiency in our case. From 
our results given in Table 2, we found no support for Hatano's result. Instead, 
the power law fit (see Fig. 2) 

Tg oc L^-^' (16) 

is almost the same as the monovariate case. This is not surprising as Berg 
mentioned that the optimal performance for multicanonical algorithm is oc 
N{=L?') based on random walk picture. 
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L 


Flat Histogram 


Multicanonical 


Bivariate in q 




(F.H.) 


(M.C.) 




4 


27.4 ± 0.14 


35.3 ± 2.8 


75.1 ± 3.0 


8 


541 ± 9 




(1.82±0.23)xl03 


12 


(4.71±0.3)xl03 


(2.61±0.45)xl03 


(9.68±1.08)xl0=^ 


16 






(3.74±0.31)xl0^ 


24 


(2.22±0.6) X 10^ 


(1.94±0.44)xl0^ 




48 




(1.46±0.52)xl06 





Table 2: Average tunneling time obtained with two dynamics for different 
lattice sizes. The multicanonical results are obtained from Ref. 0. 



The ground-state energy and entropy of the infinite system are also es- 
timated using our method. It is straightforward to obtain the ground-state 
energy in the simulation stage. We calculate the ground-state entropy from 

S{E) = ^\nn{E). (17) 

Since n{E) can be calculated from the simulation data directly, we then 
obtain S{E) with ease. 

To compare with the results obtained in the literature, we fit our data 
using the form fi = /oo + c/L^ and get e° = -1.4007 ±0.0085, s° = 0.0709 ± 
0.006. The energy fit is plotted in Fig. 3, and the entropy fit in Fig. 4. 
Our energy estimate e° = —1.4007 ± 0.0085 is consistent with the previous 
MC estimate [|] = —1.407 ± 0.008 as well as with the transfer matrix 
result § e° = -1.4024 ± 0.0012. Our entropy estimate s° = 0.0709 ± 0.006 
is also consistent with the MC estimate 0] s° = 0.071 ± 0.007 as well as the 
transfer matrix result [H] = 0.0701 ± 0.005. For the two-dimensional Ising 



spin glass system, De Simone et al. [T^ use an exact algorithm based on the 



branch- and- cut technique to find the exact ground-states with system size 
up to 50 X 50. They obtain the extrapolated result e° = -1.4022 ± 0.0003. 
When compared with Berg's result, e° = -1.394 ± 0.007, s° = 0.081 ± 0.004, 
it seems that our method gives a more accurate estimate for ground-state 
energy and entropy for an infinite system. 
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Figure 1: vs lattice size on a double log scale. 
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Figure 2: Tunneling time of bivariate model on a double log scale. 



4 Conclusions 



We have used a new approach to investigate the ground-state properties of 
the two-dimensional Ising spin glass. Compared with standard simulations, 
the advantage of our method is obvious. For the ergodicity problem en- 
countered in standard simulations, our method behaves as well as Berg's 
multicanonical ensemble method, while it is easier to be implemented com- 
pared with Berg's method. Our method also generalize straightforwardly to 
multi-variate models without much effort in programming and theory. 

To find a true ground-state, we roughly need a CPU time of order L^. 
It is the same with Lawler's exact algorithm [O. Up to size 50 x 50, De 



Simone's algorithm also needs a time of order . But it is not clear whether 
his algorithm can be efficiently implemented for 3D systems. However our 
method can also be easily applied to 3D spin glass system. If one is just inter- 
ested in finding the ground-states, there are also other optimized algorithms. 
Chen's learning algorithm [|I8| is fast in finding the ground-states compared 
with most algorithms, but it is not a general one. Thermodynamic quantities 
cannot be obtained with this algorithm. 

We believe that the approach we present in this paper is useful in studying 
the thermodynamics as well as ground-state properties for spin glass systems. 
It also can be applied to other models because of its generality. 
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